Section: New Results

Computational Statistical Physics

Participants : Grégoire Ferré, Giacomo Di Gesù, Thomas Hudson, Dorian Le Peutrec, Frédéric Legoll, Tony Lelièvre, Pierre Monmarché, Boris Nectoux, Julien Roussel, Mathias Rousset, Laura Silva Lopes, Gabriel Stoltz, Pierre Terrier, Pierre-André Zitt.

In [24], T. Lelièvre and G. Stoltz have given an overview of state-of-the art mathematical techniques which are useful to analyze and quantify the efficiency of the algorithms used in molecular dynamics, both for sampling thermodynamic quantities (canonical averages and free energies) and dynamical quantities (transition rates, reactive paths and transport coefficients).

Improved sampling methods

This section is devoted to recent methods which have been proposed in order to improve the sampling of the canonical distribution by modifying the Langevin or overdamped Langevin dynamics, or its discretization. Two general strategies have been pursued by the project-team along these lines: (i) constructing dynamics with better convergence rate and hence smaller statistical errors; (ii) the stabilization of discretization schemes by Metropolis procedures in order to allow for larger timesteps while maintaining acceptable rejection rates.

A first approach to obtaining better convergence rates consists in modifying the drift term in the overdamped-Langevin dynamics, in order to improve the rate of converge to equilibrium. This method was considered by T. Lelièvre with A. Duncan and G.A. Pavliotis (Imperial College) in [14]. It is shown that nonreversible dynamics always result in a smaller asymptotic variance (statistical error). The efficiency of the whole algorithm crucially depends on the time discretization, which may induce some bias (deterministic error). It is shown on some examples how to balance the two errors (bias and statistical errors) in order to obtain an efficient algorithm.

The discretization of overdamped Langevin dynamics, using schemes such as the Euler-Maruyama method, may lead to numerical methods that are unstable when the forces are non-globally Lipschitz. One way to stabilize numerical schemes is to superimpose some acceptance/rejection rule, based on a Metropolis-Hastings criterion for instance. However, rejections perturb the dynamical consistency of the resulting numerical method with the reference dynamics. G. Stoltz and M. Fathi (Toulouse) present in [15] some modifications of the standard stabilization of discretizations of overdamped Langevin dynamics by a Metropolis-Hastings procedure, which allow to either improve the strong order of the numerical method, or to reduce the bias in the estimation of transport coefficients characterizing the effective dynamical behavior of the dynamics.

The sampling properties of Langevin dynamics can be improved by considering more general non-quadratic kinetic energies. This was accomplished in [26], where G. Stoltz, with S. Redon and Z. Trstanova (Inria Grenoble), have studied the properties of Langevin dynamics with general, non-quadratic kinetic energies $U\left(p\right)$, showing in particular the ergodicity of the dynamics even when the kinetic force $\nabla U$ vanishes on open sets and proving linear response results for the variance of the process for kinetic energies which correspond to the so-called adaptively restrained particle simulations. This work has been complemented by [51], where G. Stoltz and Z. Trstanova provide accurate numerical schemes to integrate the modified Langevin dynamics with general kinetic energies, with possibly non globally Lipschitz momenta.

Adaptive methods

When direct sampling methods fail, it is worth considering importance sampling strategies, where the slowest direction is described by a reaction coordinate $\xi$, and the invariant measure is biased by (a fraction of) the free energy associated with $\xi$.

The first group of results along these lines concerns the study of adaptive biasing methods to compute free energy differences:

• The result obtained by H. Al Rachid (CERMICS) in collaboration with T. Lelièvre and R. Talhouk (Beirut) on the existence of a solution to the non linear Fokker Planck equation associated to the ABF process has been published, see [7].

• T. Lelièvre and G. Stoltz, together with G. Fort (Toulouse) and B. Jourdain (CERMICS), have studied the well-tempered metadynamics and many variants of this method in [41]. This dynamics can be seen as some extension of the so-called self-healing umbrella sampling method, with a partial biasing of the dynamics only. In particular, the authors propose a version which leads to much shorter exit times from metastable states (accelerated well-tempered metadynamics).

The project-team also works on adaptive splitting techniques, which forces the exploration in the direction of increasing values of the reaction coordinate. In [29], T. Lelièvre, together with C. Mayne, K. Schulten and I. Teo (Univ. Illinois), has reported on the calculation of the unbinding rate of the benzamidine-trypsin system using the Adaptive Multilevel Splitting algorithm. This is the first "real-life" test case for the adaptive multilevel splitting. In [11], T. Lelièvre and M. Rousset, in collaboration with C.E. Bréhier (Lyon), M. Gazeau (Créteil) and L. Goudenège (Centrale), propose a generalization of the Adaptive Multilevel Splitting method for discrete-in-time processes. It is shown how to make the estimator unbiased. Numerical experiments illustrate the performance of the method.

Coarse-graining and reduced descriptions

A fully atomistic description of physical systems leads to problems with a very large of unknowns, which raises challenges both on the simulation of the system and the interpretation of the results. Coarse-grained approaches, where complex molecular systems are described by a simplified model, offer an appealing alternative.

F. Legoll and T. Lelièvre, together with S. Olla (Dauphine), have proposed an analysis of the error introduced when deriving an effective dynamics for a stochastic process in large dimension on a few degrees of freedom using a projection approach à la Zwanzig [48]. More precisely, a pathwise error estimate is obtained, which is an improvement compared to a previous result by F. Legoll and T. Lelièvre where only the marginal in times were considered.

Another line of research concerns dissipative particle dynamics, where a complex molecule is replaced by an effective mesoparticle. The work [17] by G. Stoltz, together with A.-A. Homman and J.-B. Maillet (CEA), on new parallelizable numerical schemes for the integration of Dissipative Particle Dynamics with Energy conservation, has been published. Together with G. Faure and J.-B. Maillet, G. Stoltz has proposed in [16] a new formulation of smoothed dissipative particle dynamics, which can be seen as some meshless discretization of the Navier–Stokes equation perturbed by some random forcing arising from finite size effects of the underlying mesoparticles. The reformulation, in terms of internal energies rather than internal entropies, allows for a simpler and more efficient simulation, and also opens the way for a coupling with standard dissipative particle dynamics models.

G. Stoltz also suggested in [50] a new numerical integrator for DPDE which is more stable than all the previous integrators. The key point is to reduce the stochastic part of the dynamics to elementary one-dimensional dynamics, for which some Metropolis procedure can be used to prevent the appearance of negative energies at the origin of the instability of the numerical methods.

During the post-doctoral stay of I.G. Tejada (ENPC), G. Stoltz, F. Legoll and E. Cancès studied in collaboration with L. Brochard (ENPC) the derivation of a concurrent coupling technique to model fractures at the atomistic level by combining a reactive potential with a reduced harmonic approximation. The results have appeared in [28].

G. Stoltz and P. Terrier, in a joint work with M. Athènes, T. Jourdan (CEA) and G. Adjanor (EDF), have presented a coupling algorithm for cluster dynamics [52]. Rate equation cluster dynamics (RECD) is a mean field technique where only defect concentrations are considered. It consists in solving a large set of ODEs (one equation per cluster type) governing the evolution of the concentrations. Since clusters might contain up to million of atoms or defects, the number of equations becomes very large. Therefore solving such a system of ODEs becomes computationally prohibitive as the cluster sizes increase. Efficient deterministic simulations propose an approximation of the equations for large clusters by a single Fokker-Planck equation. The proposed coupling algorithm is based on a splitting of the dynamics and combines deterministic and stochastic approaches. In addition, F. Legoll and G. Stoltz have proposed in [19], with T. Jourdan (CEA) and L. Monasse (CERMICS), a new method for numerically integrating the Fokker–Planck approximation of large cluster dynamics.

Eyring–Kramers formula and quasi-stationary distributions

G. Di Gesù, T. Lelièvre and B. Nectoux, together with D. Le Peutrec, have explored the interest of using the quasi-stationary distribution approach in order to justify kinetic Monte Carlo models, and more precisely their parameterizations using the Eyring-Kramers formulas, which provide a simple rule to compute transition rates from one state to another [13]. The paper is essentially a summary of the results which have been obtained during the first two years of the PhD of B. Nectoux. A preprint with detailed proofs of these results is in preparation.

In [33], G. Di Gesù has studied with N. Berglund (Orléans) and H. Weber (Warwick) the spectral Galerkin approximations of an Allen-Cahn equation over the two-dimensional torus perturbed by weak space-time white noise. They show sharp upper and lower bounds on the transition times from a neighborhood of the stable configuration -1 to the stable configuration 1 in the small noise regime. These estimates are uniform in the discretization parameter, suggesting an Eyring-Kramers formula for the limiting renormalized stochastic PDE.

Functional inequalities and theoretical aspects

The interplay between probability theory and analysis in statistical physics is best exemplified by the functional analysis study of the semigroups associated with the generator of the stochastic processes under consideration. These generators are elliptic or hyperbolic operators. Several functional-analytic results were obtained by the team on problems of statistical physics.

D. Le Peutrec has derived Brascamp-Lieb type inequalities for general differential forms on compact Riemannian manifolds with boundary from the supersymmetry of the semiclassical Witten Laplacian [47]. These results imply the usual Brascamp-Lieb inequality and its generalization to compact Riemannian manifolds without boundary.

T. Hudson has considered with C. Hall (Oxford) and P. van Meurs (Univ. Kanazawa, Japan) the minimization of the potential energy of $N$ particles mutually interacting under a repulsive interaction potential with a certain algebraic decay assumption [42]. A major novelty of the approach is that it does not assume a finite range of interaction. The main focus of the work is on characterizing the boundary behavior of minimizers in the limit where the number of particles $N$ tends to infinity with a constant density of particles per unit volume.

G. Di Gesù has studied with M. Mariani (Rome) the small temperature limit of the Fisher information of a given probability measure with respect to the canonical measure with density proportional to $exp(-\beta V)$ [39]. The expansion reveals a hierarchy of multiple scales reflecting the metastable behavior of the underlying overdamped Langevin dynamics: distinct scales emerge and become relevant depending on whether one considers probability measures concentrated on local minima of $V$, probability measures concentrated on critical points of $V$, or generic probability measures on ${ℝ}^d$.